Monte Carlo Neutrino Oscillations 



in 
o 
o 

(N 

di- 
D 

0\ 



> : 
\o . 
in ■ 
m ■ 

o\ : 
o . 

in ■ 

o ; 
^ : 

Oh! 

i 

Oh: 

^ ■ 



X 



James P. KnelleiQ 

Department of Physics, North Carolina State University, Raleigh, North Carolina 27695-8202 and 
School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455 



Gail C. McLaughlii{[| 

Department of Physics, North Carolina State University, 

(Dated: February 7, 2008) 



North Carolina 27695-8202 



We demonstrate that the effects of matter upon neutrino propagation may be recast as the 
scattering of the initial neutrino wavefunction. Exchanging the differential, Schrodinger equation for 
an integrai equation for the scattering matrix S permits a Monte Carlo method for the computation 
of S that removes many of the numerical difficulties associated with direct integration techniques. 



PACS numbers: 02.70.Uu, 03.65.Nk, 14.60.Pq 



I. INTRODUCTION 

As a neutrino propagates through matter the non-zero 
density modulates the flavor oscillations of the neutrino 
wavefunction. The evolution of the wavefunction diffcrs 
from that in vacuum with the consequnce that neutrino 
flavour transformation may be enhanced. Appreciation 
of this effect, first discussed by Mikheyev & Smirnov and 
Wolfenstein 0, 0, is particularly important when the 
source of neutrinos is buried deep within dense matter 
such as those one can fìnd in astrophysical settings. In- 
deed, this transformation was first invoked to resolve the 
discrepancy between the observed and predicted detec- 
tion rates of solar electron neutrinos 0, Q and the most 
compelling experimental evidence for this effect has come 
from the Sudbury Neutrino Observatory which is capable 
of measuring the /iorr flavor content of the neutrinos ini- 
tially produced in the center of the Sun as electron type 
In the same fashion, the flavor content of neutrinos 
emitted from the neutrinosphere in a proto-neutron star 
will be altered by their propa gati on thro ugh the over- 
lying progenitor material 0, 15 lidi ITU IT^ . l25| and it is 
also apparent that understanding the matter effect of the 
Earth is cruciai for interpreting any future long baseline 
experiment, see e.g. 0. 

For each of these situations one is provided with the 
initial state of the neutrino at the source and wishes to 
determine the flavor content of the wavefunction after it 
has passed through the intervening material. Gauging 
the matter effects means possessing a suitable calcula- 
tional tool. The most obvious point of departure for such 
a calculation is the Schrodinger equation. Since there are 
three neutrino flavors the neutrino wavefunction must 
posses three complex components and the Hamiltonian, 
H , is a 3 x 3 matrix. In vacuum the Hamiltonian in the 
flavor basis is not diagonal and it is the presence of the 
off-diagonal terms in H that lead to flavor oscillations. 



The vacuum Hamiltonian may be diagonalized by a suit- 
able unitary transformation and it is this new basis that 
form the 'mass eigenstates'. 

In the presence of matter a potential, V[x), that takes 
into account coherent forward scattering of the neutrinos, 
must be included in the Hamiltonian. For the case of only 
active neutrino flavors (i.e. ali the flavors that have ordi- 
nary weak interactions) passing through normal matter 
the only relevant portion of V(x) is the v e — v e compo- 
nent of V(x). This is the well-known V(x) = \/2Gfp(x) 
where Gf is Fermi's Constant and p(x) is the electron 
number density. With the addition of V{x) the Hamilto- 
nian is no longer diagonal in the mass basis. A new basis, 
the 'matter eigenstates', diagonalizes H(x) but the spa- 
tial variance now within the Hamiltonian means that the 
unitary transformtion that relates the flavor to the mat- 
ter basis also varies with the propagtion distance. The 
gradient of this unitary transformation is non-zero and 
one finds that the Schrodinger equation in this new ba- 
sis picks up off-diagonal terms. Again, the presence of 
off-diagonal terms in the Schrodinger equation leads to 
mixing of the complex coefncients describing the wave- 
function and this will occur even in the matter basis if 
those terms are sufflciently large. 

Though, in general, the three complex components of 
the wavefunction oscillate simulataneously the large dif- 
ference in vacuum mass splittings usually permits us to 
consider the evolution of the neutrino wavefunction as 
being factored into two, localized, spatially separated, 
two- neutrino mixings. This factorization simplifies mat- 
ters greatly. For two-neutrino mixing there is a single 
rotation angle 6y that describes the relationship between 
the two mass eigenstates and the flavor states, and, simi- 
larly, within matter there is only one rotation angle 9(x), 
the matter mixing angle, for the relationship between 
the matter and flavor bases. In the matter eigenstate 
basis the Schrodinger equation for the evolution of the 
2-component neutrino wavefunction is 



'Electronic address: Kneller@physics.umn.edu 
^Electronic address: Gail'McLaughlin@ncsu.edu 



dx 



(IL 



k i9' 
-i9' -k 



ai 



(1) 



The prime denotes differentiation with respect to position 
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x, the quantity k is 

fc(a:) 



sin 26v 
sin 29 (x) 



(2) 



with fcy = óm 2 /(4:E) where Sm 2 is the mass squared 
difference between the neutrino mass eigenstates, and 



tan 29 (x) 



sin 26» 



v 



cos 29 v - V(x)/(2k v ) 



(3) 



Examination of equation (J2J reveals that fc passes 
through minima whenever 9 = 7r/4 and these points 
in the profile are known as the 'resonances'. The ra- 
tio 7 = k/\d9/dx\ defining the adibaticity parameter is a 
measure of the strength of the mixing and the positions 
where 7 reaches minima are 'points of maximal violation 
of adiabaticity'. It is here that the off-diagonal terms in 
equation Q are most important and the mixing is at is 
strongest. Note, as pointed out by Friedland [lj], that in 
general the positions of 'resonaces' and 'points of max- 
imal violation of adiabaticity' do not coincide and it is 
actually the latter that are more important for the evo- 
lution of the wavefunction. That said, throughout the 
remainder of the paper we will use these terms inter- 
changeably. 

The Schrodinger equation forms a starting point from 
which the neutrino wavefunction emerging from the den- 
sity profile can be determined. For 2-flavor mixing this 
is completley specified by calculating the 'survival prob- 
ability' i.e. the probability that a neutrino born as a 
particular flavor will emerge from the density profile as 
that same flavor. Unitarity provides the probability of 
detecting the counterpart flavor. Quite generally, if the 
rotation angle 9 at the neutrino source is 9 = 9q and 
the neutrino propagates to the vacuum then, after drop- 
ping the phase dependent terms, the flavor basis survival 
probability is 01 

P(v a v a ) = - [ 1 + cos 29 v cos 20 o (1 - 2P C )] . (4) 

Here the quantity Pc is known as the crossing probabil- 
ity. The crossing probability is a quantity defined in the 
matter basis and is the chance that an initial neutrino 
wavefunction transits from one matter eigenstate to the 
other. One obvious method to calculate Pc is to simply 
integrate the Schrodinger equation in the matter basis. 
If 7 is always large as the neutrino propagates then the 
off-diagonal terms in equation may be neglected, the 
integration of the Schrodinger equation is trivial and the 
wavefunction is said to evolve adibatically. There are 
also a handful of profiles where Pc has an exact ana- 
lytic solution |13l fLU , the most well-known being the 
Landau-Zener result for the infinite linear profile: 



P c = exp[-7T7c/2], 



(5) 



where j c is the adibaticity parameter evaluated at the 
resonance. The Landau-Zener equation for Pc possesses 



'troublesome pathologies' as discussed, and corrected, by 
Haxton [Ì1| . 

But exact results are scant and often one finds that nu- 
meric integration of the Schrodinger equation for many 
interesting applications can be a frustrating exercise. As 
we mentioned previously, off-diagonal terms lead to os- 
cillations and this is true even in the matter basis if the 
9' term in equation becomes large. Oscillatory So- 
lutions of differential equations obtained numerically are 
notorious for a graduai accumulation of error in both the 
phase and amplitude of the solution. A suitable change 
of variables can help to control some of these problems 
\TÌ\ but even so, with a conventional solver, one usually 
has to be very aggressive with the error control in order 
to keep the solution accurate. This requirement can lead 
to very long run times. 

In addition, the numeric integration of equation Q is 
inefficient. With a convential differential equation solver 
the increments of the integration variable (here it is x) 
are necessarily smaller than the locai oscillation length 
(~ 1/fc). In regions of very high density, such as those 
found at the centers of supernova progenitor profiles, k is 
very large and so the oscillation lengthscale will be very 
small. The differential equation solver will expend a great 
deal of time computing the wavefunction in such regions 
even though the large effective mass splitting indicates 
that the wavefunction is far from any resonance and 7 is 
large so that, in some regard, its evolution is both simple 
and uninteresting. Specialized methods, such as that by 
Petzold 9] , for highly oscillatory solutions of differential 
equations can help with the problem of small step sizes 
but their use may be limited by the requirement that any 
solution evolve adibatically. 

Due to these numeric problems, and motivated by a 
desire to comprehend the MSW effect, a number of al- 
ternate methods have been developed for calculating Pc- 
For example, one could estimate Pc by using one of the 
exact results, most typically the Landau-Zener, if that 
approximation for the profile is adequate for the situa- 
tion at hand. An alterantive approach would be to use 
the semi-analytic method by Balantekin & Beacom 0] 
for arbitrary monotonie profiles. But for one reason or 
another these alternate approaches can break down or 
are difficult to automate. One such occurence is the 
case of multiple resonances and the computational meth- 
ods used for the case of fluctuations in the solar pro- 
file m m m m are much more sophisticated than 
brute-force application of a convential differential equa- 
tion solver. 

In this paper we outline a new computational method 
for determining the neutrino wavefunction after its pas- 
sage through a density profile. Our method undertakes 
the Monte Carlo integration of a scattering matrix and 
makes no assumptions with regard to adiabaticity or 
number of resonances and devotes the bulk of the compu- 
tational time to the region around the point of maximal 
violation of adiabaticity. We derive our equations in sec- 
tion fjn] and discuss the Monte Carlo integration in sec- 
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tion £11111 We then consider the most important numerical 
difficulty in section §V\ before ending with applications 
of the technique to the density profile of the Sun and a 
density profile obtained from the evolution of supernova 
progenitor profile with a hydrodynamical calculation in 
section ilVII Throughout this paper we will only consider 
two-flavor oscillations. In an appendix we expand on our 
ponderments of practical implementations. 



II. FROM THE SCHRODINGER EQUATION TO 
THE SCATTERING MATRIX 

One persepctive on the evolution of the neutrino wave- 
function through a density profilc would be to regard the 
initial wavefunction as having been 'scattered' so as to 
produce the emerging wavefunction. The scattering ma- 
trix that relates the outgoing wavefunction to the initial 
may be derived from the Schrodinger equation. As a first 
step in its determination we define a new variable (p via 



dx 



k 

j 

7T 



so that 



= — / k dx 

7T 



(6) 
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and, secondly, we introduce a new basis, b, related to the 
matter eigenstates via 



b L 







oh 
«l 



Equation © allows us to change the independent vari- 
able from x to <f> and so measure distances in terms of 
this quantity. We see that 4> has a physical interpreta- 
tion as the number of half periods of the purely adiabatic 
solutions (Le. 9' = 0) of the Schrodinger equation QJ. 
Substitution of these definitions into equation JQ) pro- 
duces 



b H 
b L 
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b L 
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bH 
b L 



(9) 



where T is the related to the adiabaticity parameter, 



r = ir/"/. 



(10) 



Note that this definition indicates that a wavefunction 
will evolve non-adiabaticly if T >> 1. The change in 
basis allows us to focus the problem on the non-adiabatic 
part of the solution by which we mean that portion that 
jumps from one matter eigenstate to the other because, 
in this new basis, the Schrodinger equation is now purely 
off-diagonal. By integrating equation we obtain 



b H 
b L 



bm 
bw 



b H 
b L 



(11) 



(8) and repeated substitution of this result into itself yields 




where the subscripts on the -ff's indicate their argument 
with respect to </> by which we mean Hi = H(4n). This 
equation defines the scattering matrix £($) in this basis 
as 




The upper limits on the integrals appearing in equation 
(|13fl indicate the space/time ordering (0 is a monotoni- 
cally increasing function of x) but we may change ali the 



upper limits to $ by using such identities as 

/ dfa Fi / d(f> 2 H 2 = - # a / d<p 2 
{Ht H 2 0(fa - 0a) + H 2 HMfo - fa)} . (15) 

where 0(</>i — 4> 2 ) is the step function. This, and simi- 
lar identities for the the higher order multiple integrals, 
allows us to write 5*($) as 
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= 1 + (-t) 
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dt/>i / d4> 2 T{H\ H 2 ) 



3! 



(16) 



where T is the </>-ordered product. With the scattering 
matrix defined we describe our approach to its calcula- 
tion. 



III. MONTE CARLO CALCULATIONS FOR 
THE SCATTERING MATRIX 

The conversion from a differential to an integrai equa- 
tion means that completely different numerical algo- 
rithms must be applied. The number of tcrms that one 
may have to include in equation (|16fl to achieve sufR- 
cient accuracy, and the fact that H((f>) involves an e 2m ^ 
oscillatory terms, likely precludes any approach other 
than a Monte Carlo integration. Though Monte Carlo 
methods are often regarded as a last resort their use- 
fulness becomes apparent when either the boundaries of 
the integration region are very complicated or, as in this 
case, when the dimensionality of the integration measure 
means that more sophisticated algorithms will not pro- 
duce a result in a respectable amount of time. 

The quantities we select randomly are <pi to be drawn 
from a probability distribution P(4>). Naively we could 
pick values for <f>\, 4>2, ■■■ from a uniform range between 
and $ but the structure of equation (fTf)f shows that 
this would be inefhcient because the Hamiltonians are 
ali proportional to T and this quantity is largest in the 
region close to the resonance. This would suggest that 
we should select P(4>) oc T and hence use importance 
sampling for the cf>'s. This would be fine for the case of 
only one resonance but if there are multiple resonances 
we encounter problems due to the fact that T oc 8' . lì 9' 
ever switches sign then P((f>) would switch sign and, over 
some portion of the profile, we would have a negative 
probability distribution. So instead we use P(<fi) oc |T|. 

To illustrate just how sharply peaked |r|(0) can be 
we show in figure QJ this function for the BS2005- 
AGS,OP Standard Solar Model for two different 



values of sin 2 28y ■ For the upper panel we selected 
Sm 2 = 3xlCT 5 eV 2 , E = 10 Me V and sin 2 (2(9y) = 0.001 
which, as the figure indicates since the peak is |T| > 1, 
means that the resonance is non-adiabatic. The point 
of maximal violation of adiabaticity is where 7 oc 1/r 
reaches it's minimum value so by using P(4>) oc |T| as 
the probability distribution for <j> we concentrate our ef- 
forts around this point. The bottom panel shows |T| for 
the case of sin 2 29y = 0.1. For this larger value, |T| 
is less sharply peaked, the wavefunction evolves adibati- 
cally and the values of <j) we obtain from this probability 
distribution are spread over a broad range. 

Before we proceed the probability distribution must be 
normalized. The normalization A for the distribution, 
PO) = A|r|, is simply 

i/a = j #'|r| = j dx'\—\. (17) 

With P{4>) identified we can pulì out from H(<f>) the prob- 
ability distribution P(<j>) and define a reduced Hamilto- 
nian h(<p) as H{<f>) — P{4>) h(<f>); written explicitly h(<f>) 
is 

H<t>) = j {- te -**+ )■ (18) 

The definition for the scattering matrix in equation 
Ijlfijl is a sum of multiple integrals but by utilizing the 
identity 




the sum can be collapsed down to a single multiple inte- 
grai albeit one with infinite dimensionality for its mea- 
sure: 



S = (lì l iW dftj jl + (-*) hi 



+ h 2 ) + T(fti h 2 h 3 ) + ...}. ( 20 ) 



This expression for the scattering matrix is more useful tation value of the quantity s where 
from a practical standpoint because it allows us to reuse 

values of (fi. The scattering matrix is therefore the expec- -, , / \ 1 , ( — l ) nn/i ? \ 

Y 6 F s = 1 + (-«) fti + — T(ft x /i 2 ) 

+k^-T(h 1 h 2 h 3 ) + ... (21) 
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Our algorithm for the calculation of S is to simply form 
a set with samples of the matrix s and average them. 
IV. THE CROSSING PROBABILITY FROM 
THE SCATTERING MATRIX 



The scattering matrix possesses a simple structure de- 
fined by two complex numbers a and (3 so that 



S 



a (3 
-13* a* 



(22) 



It is tcmpting to regard a and (3 as Cayley-Klein pa- 
rameters but, as we shall discuss below, the Monte 
Carlo does not guarantee that S' S = 1 or, equivalently, 



a 



1 



Once S is calculated we apply S to the initial wavefunc- 
tionò(O) inorderto determine ò($), i.e&(<&) = S($)b(0). 
The crossing probabili ty, Pq, is the chance that a wave- 
function prepared in a pure matter eigenstate has tran- 
sited to the other matter eigenstate as it emerges from 
the profile. Our scattering matrix is defined in the b ba- 
sis, not the matter, a, basis and these are related by the 
expression in equation JHJ. The crossing probability is 
thus 



01 )( C e«* 



5 



(23) 
(24) 



FIG f : The absolute value of the inverse adiabaticity param- 
eter F as a function of the period counting variable <f>. The 
density profile is the BS2005-AGS,OP Standard Solar Model 
[23| and we selected Sm 2 = 3 x fO" 5 eV 2 , E = 10 MeV. In 
the top panel we used s\v?{2 8v) = 0.001 and for the bottoni 
sin 2 (26V) = 0.1. 



The superscript upon is to remind the reader of the 
second line of this equation. 

We may also define Pc as being the difference from 
unity of the probability that a wavefunction prepared in 
a pure matter eigenstate survives as that same matter 
eigenstate as it emerges from the profile: i.e. 



Pg = 1 - ( 1 ) 5 f 







(f 0) 



-tir® 







s 



(25) 
(26) 
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Again, the superscript upon Pg is to remind the reader 
of the second line of this equation. 

Thus our scattering matrix can be used to construct 
two values for Pc and if S were unitary then they would 
be equal. The difference between them is due to finite 
sampling and is the subject of the next section. 



V. THE DISTRIBUTIONS OF Pg AND Pg AND 
THE UNITARITY OF S 

After execution of the Monte Carlo algorithm for a 
given profile and mixing parameters, one obtains a scat- 
tering matrix S from which Pg and Pq can be formed. 
The scattering matrix calculated by this method does 
not, in general, guarantee that the identity S — 1 = 
is satisfied. This is equivalent to the statements that 
\a\ 2 + |/3| 2 - 1 f and - Pg f 0. Thus a and 
are not Cayley-Klein parameters and the two calculated 
crossing probabilities are not exactly equal. Also, if the 
calculation for a given profile and mixing paramaters is 
repeated then we obtain a new scattering matrix and new 
values for Pg and Pq. The difference between the two 
crossing probabilities for a given run and their change 
from one run to the next arises because we only con- 
struct a finite set of samples of s. Only in the limit of an 
infinite number of samples would Pg and Pq be exactly 
equal and our calcualtion give the same result every time. 

We stress that this behavior is not a fundamental fìaw 
of the Monte Carlo technique but rather a numeric issue 
related to the usuai lack of infinite computing resources. 
For this reason one must be content with values for Pg 
and Pq that differ from the true crossing probability and 
from each other. With the cautionary note that what 
follows is specific to our implementation of the algorithm 
and the test problems we selected, we try and provide 
some guidance on how to obtain the most accurate cal- 
culation in the least computational timo. 

The values of Pg and Pq obtained from a given calcula- 
tion are drawn from parent distributions that, in general, 
are unique to the particular profile, mixing parameters 
and also the implementation of the algorithm. Thcse par- 
ent distributions may be reconstructed by repeating the 
calculation for Pg and Pq until a sufficiently large sam- 
ple of results has been extracted. As an example, the fre- 
quency distributions of Pg and Pq for the case of 10 MeV 
neutrino passing through the BS2005-AGS,OP Standard 
Solar Model density profile with Sm 2 — 3 x 10 " 5 eV 2 
and sin 2 (2 6V) = 0.001 are shown in figure The ini- 
tial impression is that the distributions for Pg and Pq 
are both Gaussian with similar means and variances and 
indeed the <i-statistic from a Kolmogorov-Smirnov test 
(using the Lilliefors 15f criticai values) verifies this con- 
clusion. The departure from unitarity is affected by the 
number of trials that go into the calculation of S. In 
figure © we plot the mean value of S't S — 1 and the 
1 — a spread of the sample for various values of Nt- 
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FIG. 2: The frequency distribution of Pg (solid) and P§ 
(dashed) of 10,000 results from the Monte Carlo calcula- 
tion using the BS2005-AGS,OP Standard Solar Model |H, 
Sm 2 = 3 x IO" 5 eV 2 , E = 10 MeV and sin 2 (2 6V) = 0.001 as 
the physical parameters. The number of trials is Nt = IO 4 . 



Again, the calculation is for a neutrino passing through 
the BS2005-AGS,OP Standard Solar Model [23| density 
profile and for Sm 2 = 3 x IO" 5 eV 2 , E = 10 Me V and 
sin 2 (2 6» v ) = 0.001. The mean values of SÌ S - 1 ali lie 
above zero indicating that the mean value of Pq is appar- 
ently slightly larger than the mean for Pg but that this 
difference disappears as Nt increases. The error bars on 
each point are the 1 — a spread in S< S — 1 and these 
clearly diminish as Nt increases. We find that the size 
of the error bars follows a trend proportional to 1/^/Nt 
as indicated by the dashed lines in the figure. The figure 
indicates unitarity is achieved in the limit when the num- 
ber of samples of the scattering matrix becomes infinite. 

But Gaussianity in the distributions for Pg and Pq 
does not always occur and should not be taken for 
granted. If we consider a case where Pc is dose to 
zero non-Gaussiani ty becomes apparent. In figure |0J we 
show the frequency distributions for the case of Sm 2 = 
3xl0" 5 eV 2 ,£= 10 MeV and sin 2 (2 8 V ) = 0.1. The dis- 
tribution for Pg remains Gaussian but we immediately 
notice that the distribution for Pq has changed and we 
find that a Gamma distribution with an a parameter that 
is close to unity is a good fit. The figure also shows that 
negative values Pg can be obtained whereas Pq is always 
positive. This result arises due to the definitions in equa- 
tion (|24|l and (|26[) : values of Pq less than zero are not 
allowed but there is no similar restriction for Pg. Note 
also the very different widths of the distributions: the 
width of the Pg distribution is similar to that in figure 

@ but the width of Pq has shrunk considerably. This 
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FIG. 3: The mean value of S — 1 indicateci by the squares 
as a function of Nt for the calculation outlined in the text. 
The error bars are not the error in the mean but rather indi- 
cate the 1 — a spr ead in the values of Sft S — 1. The dashed 
lines are St S -1±C/VNt with C = 0.468. 



FIG. 4: The frequency distribution of Pg (solid) and Pg 
(dashed) of 10,000 results from the Monte Carlo calcula- 
tion using the BS2005-AGS,OP Standard Solar Model Q, 
Sm 2 = 3 x IO" 5 eV 2 , E = 10 MeV and sin 2 (2ÉV) = 0.1 
as the physical parameters. The number of trials is again 
N T = IO 4 . 



difference in the widths of the two crossing probabilitics 
would indicate that the deviation from unitarity, S—l, 
will be dominated by the spread in Pg, the values of Pq 
having such a small variance. Thus, when Pq is close 
to zero Pq is much more accurately calculated than Pg. 
We also find that for this test case, the width of the two 
distributions varies with Nt in different fashions. For Pg 
the spread again varies as 1/^/Nr but the width of the 
Pq distribution now behaves as 1/Nt- From additional 
test cases we found that that when Pq approaches unity 
it is Pg that is the more accurately calcualted. Our ex- 
perience has also shown that in some circumstances the 
distributions for Pg and Pq can also change shape as 
Nt is varied: for small Nt the distibution may be like 
a Gamma distribution with a modest a parameter but 
then will morph to something closer to a Gaussian dis- 
tribution as Nt increases. 

These results hint at the interesting underlying numer- 
ics of this Monte Carlo approach but they also introduce 
some confusion into what would be a reasonable modus 
operandi. The parent distributions for Pg and Pq are 
not, in general, the same and we do not know a priori 
their shape or if they are similar. This would seem to 
preclude combining the results for Pg and Pq in some 
way so as to obtain a more accurate result. The accuracy 
of the results depend upon Nt but in a way that varies 
as we change the profile and mixing parameters. Before 
we do the calculation we do not know how large we must 
make Nt to reach our intended level of accuracy. In prac- 
tice we adopted a 'worst case scneario' approach whereby 



we calculate both Pg and Pg, assuming that the accuracy 
varies as 1/ ' \[N~t. One would then require Nt ~ IO 6 tri- 
als to reach a level of accuracy of ~ 0.1%. We then used 
Pq for the crossing probability if Pc < 0.5 and Pg oth- 
erwise. As we said, the shape of the distributions for Pg 
and Pq can vary with the number of samples so breaking 
up the Nt trials into a number of smaller runs (e.g. 10 
runs with IO 5 samples in each), calculating Pg and Pq 
from each run and then averaging the results must also 
be approached with caution. To avoid potential bias in 
such a procedure we only accepted the result from the one 
run with the full number of trials we specified. Though 
this conservative approach has the drawback that the 
runtime of our code may be longer than necessary the 
results always achieve our desired level of accuracy and 
often considerably so. 

VI. EXAMPLE CALCULATIONS 

We finish with three applications of our method. We 
first demonstrate the method with a calculation of the 
survival probability of electron neutrinos using the solar 
density profile and two different values of sin 2 26>y . We 
then go on to use profile from an aspherical supernova 
simulation, which involves multiple resonances. 

The passage of neutrinos through the solar density pro- 
file is a well studied problem and therefore there are 
a number of already published calculations. In figure 
we calculate the survival probability of electron neu- 
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FIG. 5: The neutrino survival probability through the Stan- 
dard Solar Model density profile as a function of 5m 2 /E. In 
the top panel sin 2 (2#y) = 0.001, in the bottoni sin 2 (28v) ~ 
0.1. The source of neutrinos is located at 3/10 Rq and the 
neutrinos propagate back through the core and out the other 
side. The error bars on each point are the rms spread in the 
results from 8 repetitions. The Gaussian estimator leads to a 
bias so the accuracy should be regarded as illustrative. 



trinos over a spectrum in energy through the BS2005- 
AGS,OP Standard Solar Model [23J. For these figures 
we select either sin 2 (2 6> y ) = 0.001 or sin 2 (2 6V) =0.1. 
The source of neutrinos is located at 3/10 of the solar 
radius and they propagate back through the core and 
emerge the other side. These figures agree those of Hax- 
ton [16( for the sanie calculation. In these calculations 
the lower energy neutrinos experience a doublé resonance 
while the higher energy neutrinos experience only one. 
This changeover is seen in the bottoni panel of figure 

where the survival probability transits from ~ to 
1 at Sui 2 /E ~ 10~ 6 eV 2 /MeV. The top panel in figure 

exhibits rapid fìuctuations in the survival probabil- 
ity (which are by no means resolved with energy spacing 
we used) and indicate phase effects as discussed in fl7| . 
These features in the figure demonstrate that the Monte 
Carlo is capable of reproducing the results of other cal- 
culations. 



The measured value of dy, solar is larger than what we 
have used in our example calculations here, and there- 
fore neutrinos from the sun go through adiabatic neu- 
trino flavor transformation. However, the value of #13 
is yet unknown. This angle will determine the degree 
of flavor transformation in the core collapse supernova. 

Therefore, we consider finally the more complicated 
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FIG. 6: The electron neutrino potential energy, V = 
function of radiai distance for the model dis- 
cussed in the text. The upper figure is the entire profile, the 
lower focuses upon that portion up to 6 x IO 9 cm. In both pan- 
els the dashed lines indicate the resonance potential energies 
for 5.4 MeV (upper) and 16 MeV (lower) neutrinos indicating 
that neutrinos with energies between these values will expe- 
rience a triple resonance. The mass splitting is chosen to be 



Sm 2 



3 x 1Q- S eV 2 and sin 2 (26 v ) = 4 x 10" 



profile shown in figure (J5J). This profile is a product 
of the evolution of a supernova progenitor model using 
the VHI hydrodynamical code. An l = 2 spherical har- 
monic velocity perturbation was inserted by hand into 
the ul3.2 progenitor profile of Heger [24| to cause the 
star to explode asymmetrically. As a consequence of the 
asphericity several density minima were produced and 
the profile shown is a radiai slice through the model 9 s 
after the bounce. We select 8m 2 — 3 x IO -3 eV 2 and 
sin 2 (2(V) = 4 x 10~ 4 and find that neutrinos with ener- 
gies between 5.4 MeV and 16 MeV will experience a triple 
resonance, this region is magnified in the lower panel of 
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FIG. 7: The neutrino survival probability through the density 
profile shown in figure © as a function of the neutrino energy. 
For this calculation Sm 2 = 3 x 1CT 3 eV 2 and sin 2 (2#v) 
'. x 10~ 4 . Again, the error bars on each point are the rms 
spread in the results from 8 repetitions of the calculation using 
the Gaussian estimator. Though this estimator is biased they 
indicate the accuracry of the result. 



figure ijHJ- The results of the calculation are shown in 
figure (JZJ). Between 5.4 MeV and 16 MeV the survival 
probability, again, exhibits phase effeets. At present we 
only wish to illustrate a potential use of the technique 
presented here for the case of multiple resonances. The 
calculation leading to this profile will be discussed else- 
where along with a more detailed studied of the observ- 
able consequences of multiple resonances from aspherical 
supernova explosions [25| . 



VII. SUMMARY 

We have shown that the effeets of matter upon the 
propagation of neutrinos may be described as the scat- 
tering of an initial neutrino wavefunction permitting us 
to exchange the differential Schrodinger equation for an 
integrai equation for the scattering matrix. In this for- 
mulism we are able to avoid the numerical difficultics 
associated with oscillatory solutions to differential equa- 
tions by, instead, using Monte Carlo integrators that fo- 
cus the calculation onto the most important aspeets of 
the problem. Though slow to converge compared to more 
sophisticated methods, and possessing inherent numeri- 
cai error due to finite sampling, Monte Carlo integrators 
have the advantage of easily controlied runtimes and the 
numerical errors are both understood and ameliorated by 
repctition. This technique may be useful in a number of 
interesting density profiles that are difficult to work with 
using traditional techniques, particularily those involv- 
ing multiple resonance regions, such as in a core collapse 
supernovae. 
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APPENDIX A: PRACTICAL CONSIDERATIONS 

As we described in section flllll the scattering matrix 
is found to be the average of a set of samples for s. As a 
reminder, s is given by 

s = 1 + (-») fe x + T(/n fea) 

+ i^t T (h 1 h 2 h 3 ) + ... (Al) 

OC 

i=0 

and the subscripts on the reduced Hamiltonians h indi- 
cate the (f> argument by which we mean hi = h(4>i) and 
h(<p) is given by equation (|18J) . Constructing a set of s to 
average is the principle task of the algorithm. It is not 
our intention to proscribe a recipe for the construction of 
s, and the reader can find many additional runtime sav- 
ings that are not discussed here, but rather we outlinc 
some considerations that may be useful. 

1. Truncating the series 

Formally s is the sum of an infinite number of terms 
but in practice we must truncate the series at some order 
Ns- The basis for selection of Ns comes from noticing 
that the terms in s are proportional to a unitary matrix 
and a weighting factor w with 



We can set a value for Ns by requiring that the weight of 
the terms we retain are larger than some specified level 
£5; that is, wn s > es- Since the weights are inversely 
related to the normalization the smaller the value of A 
then then larger Ns- Small value of A, as seen in equa- 
tion H17JI . occur for greater differences between the initial 
and final rotation angles across a resonance and/or the 
greater the number of zeros for T. The value of es should 
be sufficiently small that the numerical error in the values 
of Pq and P^, should be dominated by the finite sampling 
error otherwise the crossing probabilities would contain 
a systematic error due to this truncation. For a desired 
level of accuracy of - 0.1% in Pg and P§ we found that 
£5 — 10~ 4 was sufficient. 

2. Generating random values of <j> 

With Ng chosen the structure of equation i|Al|l indi- 
cates we need Ns values of <f> to compute s. The most 
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efficient method for obtaining a sequence of <j>'s from the 
probability distribution P(4>) is to relate P((f>) to the uni- 
form distribution so that one may use a pseudo-random 
numbcr generator. This is achieved by calculating the 
accumulated probability, F{4>), from 

F(4>)= [ A\T{<f/)\d4f (A4) 

and then inversion of the relationship to form 4>{F). The 
requirement that F($>) — 1 sets the normalization A as 
shown in equation ljl7Jl . After substituting the definition 
of T, equation (TDJ, and </>, equation ©, we find that for 
a monotonie profile F(<p) = A\9(4>) ~ 9(0)\. This result 
suggests that for a general profile we can also avoid per- 
forming the integration if we identify the zeros of T and 
break apart the profile at those points so as to create a 
series of monotonie profiles. The absolute difference of 
the mixing angle across each monotonie region can be 
computed and the calculation for F((f>) is then an appro- 
priate summation. The advantages of calculating F(cj>) 
this way are: firstly, that it is far quicker than doing the 
integration, and secondly, T can be somewhat noisy - as 
shown in figure Q - due to numerical problems associ- 
ated with forming derivatives. 

To use the relationship between <p and F one generates 
a pseudo-random number u from a uniform probability 
distribution and sets F = u before inserting this value 
into 4*{F)- There is one circumstance where inversion 
of F((f>) to 4>{F) is not possible and this occurs when- 
ever |T| oc \9'\ = over some extended distance within 



a profile. Such a region would possess a Constant den- 
sity. But over this region S = 1 in the b basis so there 
is no need to perform the Monte Carlo calculation for 
this region. If this situation arises a simple solution is, 
again, to break apart the profile and only calculate the 
scattering matrix for those regions where |T| ^ 0. In this 
way the total scattering matrix for the entire profile is 
the ordered product of the scattering matrices for each 
|r| ^ zone. 

3. Efficiently using the random <f> 

Once the Ag values of <j> have been found and stored 
in an array, a possible algorithm for s would be: 

1. use the first value, (j>\, to calculate s\ and add it to 
the unit matrix, 

2. 0-order the first two values, <p\ and 02, calculate 
S2, and add it to the 1 + si sum, 

3. repeat for ali Ag terms. 

In this scheme each term in s is calculated just once. But 
the presence of the weighting factors Wi indicate that this 
is not optimal: we should calculate a term Si much more 
frequently if its weight is large and less frequently if the 
weight is small. There are many ways one can achieve 
a better load balancing: we adopted, after realizing that 
the labels on the 0's may be swapped amongst them- 
selves, a scheme whereby we rewrite equation (|Al(l as 



s = 1 



iV 

E 



t± h . 

Ci J 



N, 



H) 2 



E 



T(hj h k ) + ...+ 



H) 



N^ + l 



(A^ + l)! 



T{hi h 2 h 3 ... h N + i) + 



(A5) 



where at^C are the binomial coefficients, is an integer 
that satisfies 1 < < Ns and {j, k} indicates ali com- 
binations of two (f>'s from the first Nj> in the list. This 
equation expresses the fact that any element of the first 
A^ values of <fi from our array may be used to calculate 
Si, any ordered pair of the first A^ may be selected for 
S2 and so up to sa^, thereafter we calculate the higher 
order terms as described before. The appearance of the 
binomial coefficients in the denominators has the effect 
of increasing the number of trials that will form the first 
A^ terms of S. But the additional computation obvi- 
ously leads to an increase in the amount of time required 
to generate just one s. To compensate for the longer 
runtime we can reduce the number of samplcs of s that 
we average to form the scattering matrix. If t(0) is the 
time required to form s via equation (|A1I) . and r(N ( f,) 
is the amount of time to calculate s according to equa- 



tion 1|A5(I . then the number of s samples that we would 
have averaged with equation I|A1|I . which we cali At(0), 
is reduced to Nt(N,p) when we use equation HA5() so as 
maintain N T (N,p) T(A r ) = Ar(0) r(0). Even though the 
number of s that we average to form the scattering matrix 
is reduced a judicious choice for and the presence of 
the binomial coefficients can more than compensate this 
loss so that our scattering matrix is more accurate and 
the code more efficient. We base our decision for selecting 
A^ by defining a quantity Vg as 



V s = 



w 



' I n^Cì + 



N s 

ì=Na+1 



W? 



r(0) 



and determine the value of A^ that minimizes Vs- The 
reader may find that an alternate selection criteria leads 
to a more efficient algorithm. The computation times r 
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were found by numerical experiments and the applica- alternatively have some knowledge of their relative size 
tion of fitting formulae to the results although one may from the design of the algorithm. 
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